//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : CudaWavesComputer.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef CUDAWAVESCOMPUTER_H
#define CUDAWAVESCOMPUTER_H

#include "defines.h"
#include "Array1D.h"
#include "Array2D.h"
#include "CUDAPACKAGE/common/globals.h"

//----------------------------------------------------------------------------
// class CudaWavesComputer
//----------------------------------------------------------------------------
class CudaWavesComputer
{

public:

    CudaWavesComputer(const int Nx1,
                      const int Nx2,
                      const REAL dx1,
                      const REAL dx2,
                      const int bcN,                      
                      const int bcE,      
                      const int bcS,      
                      const int bcW,      
                      const REAL g);

    virtual ~CudaWavesComputer();

    void initialize(const Array2D<REAL> *hU1,
                    const Array2D<REAL> *hU2,
                    const Array2D<REAL> *hh0,
                    const Array2D<SHORT> *hwet,
                    const Array2D<REAL> *hD0,
                    const Array2D<REAL> *hSRCz,
                    const Array2D<REAL> *hSRCp,
                    const Array2D<REAL> *halfa_brk,
                    const Array2D<REAL> *hzeta_std,
                    const Array2D<REAL> *hSppC,
                    const Array1D<REAL> *hstoN,
                    const Array1D<REAL> *hstoE,
                    const Array1D<REAL> *hstoS,
                    const Array1D<REAL> *hstoW);

    void shiftZetaPhi(const REAL eps);

    void waveBreaking(Array2D<REAL> *hzetastd,
                      Array2D<REAL> *halfabrk,
                      REAL zeta_mean, REAL tau_brk, REAL dt_new);

    void boundaryZetaPhi(Array2D<REAL> *hphi,
                         Array2D<REAL> *hzeta,
                         const Array1D<REAL> *hphi_in_N_ext,
                         const Array1D<REAL> *hphi_in_N_cen,
                         const Array2D<REAL> *hphiN,
                         const Array1D<REAL> *hzeta_in_N_ext,
                         const Array1D<REAL> *hzeta_in_N_cen,
                         const Array2D<REAL> *hzetaN,
                         const Array1D<REAL> *hphi_in_E_ext,
                         const Array1D<REAL> *hphi_in_E_cen,
                         const Array2D<REAL> *hphiE,
                         const Array1D<REAL> *hzeta_in_E_ext,
                         const Array1D<REAL> *hzeta_in_E_cen,
                         const Array2D<REAL> *hzetaE,
                         const Array1D<REAL> *hphi_in_S_ext,
                         const Array1D<REAL> *hphi_in_S_cen,
                         const Array2D<REAL> *hphiS,
                         const Array1D<REAL> *hzeta_in_S_ext,
                         const Array1D<REAL> *hzeta_in_S_cen,
                         const Array2D<REAL> *hzetaS,
                         const Array1D<REAL> *hphi_in_W_ext,
                         const Array1D<REAL> *hphi_in_W_cen,
                         const Array2D<REAL> *hphiW,
                         const Array1D<REAL> *hzeta_in_W_ext,
                         const Array1D<REAL> *hzeta_in_W_cen,
                         const Array2D<REAL> *hzetaW);

    void solvePsi(Array2D<REAL> *hbpsi,
                  const REAL hpslow,
                  const Array2D<REAL> *hpsiN,
                  const Array1D<REAL> *hpsi_in_N_ext,
                  const Array1D<REAL> *hpsi_in_N_cen,
                  const Array2D<REAL> *hpsiE,
                  const Array1D<REAL> *hpsi_in_E_ext,
                  const Array1D<REAL> *hpsi_in_E_cen,
                  const Array2D<REAL> *hpsiS,
                  const Array1D<REAL> *hpsi_in_S_ext,
                  const Array1D<REAL> *hpsi_in_S_cen,
                  const Array2D<REAL> *hpsiW,
                  const Array1D<REAL> *hpsi_in_W_ext,
                  const Array1D<REAL> *hpsi_in_W_cen);

    void adjustPsi(const REAL dt_new,
                   const REAL dt_old);

    void shiftPsi();

    void boundaryPsi(Array2D<REAL> *hpsi,
                     const Array2D<REAL> *hpsiN,
                     const Array2D<REAL> *hpsiE,
                     const Array2D<REAL> *hpsiS,
                     const Array2D<REAL> *hpsiW);

    void timeDerivatives(const REAL hpslow,
                         const Array2D<REAL> *hpressurePulse);

    void smoothZetaPhi(const REAL eps);

    void leapFrog(Array2D<REAL> *hzeta,
                  const REAL dt_new,
                  const REAL dt_old);

    REAL* getDevicePtr_bpsi();

    REAL* getDevicePtr_psi();

private:

    /// @cond nevercall
    CudaWavesComputer();
    CudaWavesComputer(const CudaWavesComputer& m);
    CudaWavesComputer & operator =(const CudaWavesComputer& m);
    /// @endcond

    void setup();

    void cleanup();

    void determineSizes();

    const int m_Nx1;
    const int m_Nx2;
    const REAL m_dx1;
    const REAL m_dx2;
    const int m_bcN;
    const int m_bcE;
    const int m_bcS;
    const int m_bcW;
    const REAL m_g; // gravity

    int m_Nc;
    int m_Nr;
    int m_n;
    int m_numBlocks;
    int m_numThreads;
    size_t m_numBytes;


    // device data objects
    REAL *m_dpsi;          // psi
    REAL *m_dpsim1;        // psi 1 step earlier
    REAL *m_dpsim2;        // psi 2 steps earlier
    REAL *m_dpsi_in_N_ext;
    REAL *m_dpsi_in_N_cen;
    REAL *m_dpsiN;
    REAL *m_dpsi_in_E_ext;
    REAL *m_dpsi_in_E_cen;
    REAL *m_dpsiE;
    REAL *m_dpsi_in_S_ext;
    REAL *m_dpsi_in_S_cen;
    REAL *m_dpsiS;
    REAL *m_dpsi_in_W_ext;
    REAL *m_dpsi_in_W_cen;
    REAL *m_dpsiW;
    REAL *m_dphi_in_N_ext;
    REAL *m_dphi_in_N_cen;
    REAL *m_dphiN;
    REAL *m_dphi_in_E_ext;
    REAL *m_dphi_in_E_cen;
    REAL *m_dphiE;
    REAL *m_dphi_in_S_ext;
    REAL *m_dphi_in_S_cen;
    REAL *m_dphiS;
    REAL *m_dphi_in_W_ext;
    REAL *m_dphi_in_W_cen;
    REAL *m_dphiW;
    REAL *m_dzeta_in_N_ext;
    REAL *m_dzeta_in_N_cen;
    REAL *m_dzetaN;
    REAL *m_dzeta_in_E_ext;
    REAL *m_dzeta_in_E_cen;
    REAL *m_dzetaE;
    REAL *m_dzeta_in_S_ext;
    REAL *m_dzeta_in_S_cen;
    REAL *m_dzetaS;
    REAL *m_dzeta_in_W_ext;
    REAL *m_dzeta_in_W_cen;
    REAL *m_dzetaW;
    REAL *m_dstoN;
    REAL *m_dstoE;
    REAL *m_dstoS;
    REAL *m_dstoW;
    REAL *m_dbufc;
    REAL *m_dbufr;
    REAL *m_dpressurePulse;
    REAL *m_dU1;
    REAL *m_dU2;
    REAL *m_dh0;
    SHORT *m_dwet;
    REAL *m_dD0;
    REAL *m_dzeta;
    REAL *m_dzetam1;
    REAL *m_dzetam2;
    REAL *m_ddt_zeta; // derivative d zeta / dt
    REAL *m_ddt_zetam1;
    REAL *m_dphi;
    REAL *m_dphim1;
    REAL *m_dphim2;
    REAL *m_ddt_phi;  // derivative d phi / dt
    REAL *m_ddt_phim1;
    REAL *m_dSRCz;
    REAL *m_dSRCp;
    REAL *m_dSppC;
    REAL *m_dbpsi;
    REAL *m_dalfa_brk;
    REAL *m_dzeta_std;

    unsigned int m_nx[20];
    unsigned int m_ny[20];
    unsigned int m_cx[20];
    unsigned int m_cy[20];
};

#endif // CUDAWAVESCOMPUTER_H
